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ABSTRACT 


A  new  technique  for  estimation  of  the  spectrum  of  an 
unknown  signal  is  investigated.  The  technique  combines  a 
radiometric  method  of  signal  detection  (modulating  the  re¬ 
ceiver  noise  with  the  assumed  signal  by  alternately  connect¬ 
ing  and  disconnecting  the  input  transducer)  with  a  Fast 
Fourier  Transform  processor  for  spectral  analysis.  After 
obtaining  averaged  spectral  estimates  with  the  transducer 
(antenna)  connected  and  disconnected  the  difference  is  formed 
and  statistical  analyses  applied  to  determine  statistically 
significant  terms.  The  significant  components  comprise  the 
estimate  of  the  spectrum  of  the  signal. 

Analytical  results  are  derived  for  sinusoidal  signals  and 
uncorrelated,  Gaussian  noise.  Experimental  results  for  both 
narrow  and  wideband  signals  and  additive  Gaussian  noise  are 
presented.  Reasonable  spectral  estimates  of  a  wideband  signal 
are  obtained  with  16  averages  when  the  composite  signal-to- 
noise  ratio  is  0  dB. 
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I .  INTRODUCTION 


Wideband  noise-like  signals  occur  in  various  types  of 
radio  communications  and  radio  ranging  systems.  In  this  re¬ 
port,  we  present  a  method  of  signal  processing  which  provides 
a  measure  of  the  spectrum  of  a  wideband  signal.  This  method 
combines  the  technique  of  radiometry  which  is  used  to  detect 
weak  radio  signals  with  the  Fourier  transform  which  defines 
the  spectrum  of  an  electrical  signal. 

For  more  than  twenty  years,  radiometers  have  been  used  in 
radio  astronomy  and  elsewhere  to  detect  low  power  signals 
(weak  stars) .  In  a  radiometer  the  receiving  element  (antenna) 
is  alternately  connected  to  and  disconnected  from  the  receiver 
as  shown  in  Fig.  1.  The  received  voltage  is  therefore  alter¬ 
nately  weak  signal  s(t)  plus  noise  n(t)  and  noise  only.  Hence, 
the  receiver  noise  serves  as  a  "carrier"  which  is  then  ampli¬ 
tude  modulated  by  the  weak  signal  by  virtue  of  the  switching 
action.  Typically,  a  square-law  device  is  used  to  demodulate 
the  carrier.  The  square-law  detector  makes  the  system  a  power 
detector.  Integration  provides  energy  detection  capability. 

The  output  from  the  energy  detector  alternately  represents 
the  energy  level  with  no  signal  present  (antenna  disconnected) 
and  the  energy  level  when  it  is  possible  that  a  signal  is 
present  (antenna  connected) .  The  two  values  are  compared 
using  appropriate  statistical  analysis  and  a  decision  is 
made  concerning  the  presence  or  absence  of  a  signal.  iRef.  1] 
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So,  the  use  o£  the  receiver's  own  noise  as  a  reference 
makes  possible  the  detection  of  weak  signals  when  the  system 
noise  is  predominantly  receiver  noise  (thermal  noise-limited 
receivers) .  [In  general,  systems  are  thermal-noise  limited 
if  the  frequencies  of  interest  exceed  20  MHz.]  Our  interest 
is  in  the  measurement  of  the  spectra  of  weak  signals.  Re¬ 
placing  the  square-law  detector  of  Fig.  1  with  a  Fourier 
transform  processor  results  in  the  system  of  Fig.  2.  (In  a 
practical  system,  the  receiving  subsystem  of  Fig.  2  includes 
a  mixer  so  that  the  FFT  processor  operates  on  a  lowpass  or 
baseband  portion  of  the  spectrum.) 

In  this  report,  the  processor  of  Fig.  2  includes  a  sampler 
which  then  allows  use  of  the  Fast  Fourier  Transform  (FFT) 
algorithm  to  generate  a  discrete  form  of  the  spectrum  of  the 
processor  input  voltage.  The  FFT  processor  output  is  manipu¬ 
lated  to  produce  an  array  of  numbers  proportional  to  the  power 
at  discrete  frequencies,  (the  discrete  power  spectrum  of  the 
input) .  An  average  over  M  blocks  of  input  data  is  formed  which 
provides  an  output  proportional  to  the  energy  of  the  input. 

Because  of  the  switching  action  at  the  system  input,  the 
FFT  processor  output  alternately  represents  the  spectrum 
^(f)  of  the  noise  and  the  spectrum  G^,( f)  of  the  combination 
of  signal  and  noise.  For  the  case  of  additive  noise  uncorre¬ 
lated  with  the  signal,  we  can  write  the  signal  spectrum  Gg(f) 
as  Gc(f)  -  GN(f).  Therefore,  by  subtracting  the  averaged 
spectra  obtained  with  the  FFT  processor,  a  measure  of  the 


spectrum  of  the  signal  is  formed  when  a  signal  is  present. 

A  near-zero  difference  indicates  the  absence  of  a  signal. 

Fig.  3a  shows  the  spectrum  GN(f)  of  the  noise  only.  Fig.  3b 
shows  the  spectrum  Gc(f)  of  the  combined  signal  and  noise  for 
a  square  wave  signal,  and  Fig.  3c  shows  the  difference  Gg(f). 
Statistical  analyses  relate  the  number  of  blocks  of  data  used 
in  averaging  to  false  alarm  rate  and  detection  probabilities. 

In  this  manner,  the  significant  values  of  the  difference 
(signal)  spectrum  are  determined.  Fig.  3c  shows  that  the 
seventh  and  higher  order  harmonics  of  the  square  wave  are 
missing  from  the  spectral  estimate.  This  is  because  the 
energy  content  of  the  higher  order  harmonics  cannot  be  dis¬ 
cerned  from  the  energy  content  of  the  noise.  Fig.  3c  also 
shows  two  components  which  are  not  part  of  the  spectrum  of 
the  signal  (one  at  approximately  30  and  another  at  approxi¬ 
mately  150) .  These  are  caused  by  false  alarms  in  the  statis¬ 
tical  analysis. 

In  this  report,  we  identify  this  method  with  the  acronym 
FABLUS  (Fourier  Analysis  of  Broadband,  Low-level  Signals 
Using  Switching) . 

The  FABLUS  method  of  estimating  signal  spectra  is  useful 
when  the  receiver  noise  limits  system  performance.  Noise  which 
appears  at  the  antenna  (transducer)  terminals  is  treated  as 
signal.  For  accurate  estimates  of  signal  spectrum  this  method 
also  requires  the  noise  to  be  stationary  over  the  interval  of 
time  required  to  obtain  all  the  data  used  in  calculating  the 
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averages  of  the  spectra  of  noise  and  combination  of  signal  and 
noise.  When  these  two  conditions  are  satisfied,  there  is  no 
theoretical  limit  on  the  level  of  a  signal  which  can  be  de¬ 
tected  and  the  accuracy  with  which  the  spectrum  can  be  mea¬ 
sured.  By  averaging  over  a  large  number  of  blocks  of  data 
(many  ensembles  or  many  cycles  of  the  switch) ,  the  spectrum 
of  a  comparatively  very  weak  signal  can  be  accurately  measured 
by  this  method. 

Table  1  is  an  example  of  the  results  obtained  using  this 
method.  Shown  is  the  probability  of  detecting  a  signal  com¬ 
ponent  as  a  function  of  the  ensemble  size  M  and  the  ratio  of 
the  average  signal  amplitude  to  the  average  of  the  noise 
amplitude  (frequency  interval  signal-to-noise  ratio)  when  the 
false  alarm  rate  is  2.5  percent. 

The  concept  of  signal-to-noise  ratio  requires  definition 
when  using  the  discrete  Fourier  transform.  For  example,  if 
one  watt  of  white,  Gaussian  noise  (bandlimited  to  the  pro¬ 
cessor  bandwidth)  is  applied  to  an  FFT  processor  which  uses 
a  1024  sample  block  length  the  output  will  contain  1/512  watts 
in  each  frequency  interval.  If  a  one  watt  sinusoidal  signal, 
centered  on  one  of  the  frequency  intervals,  is  applied  to 
the  processor  the  output  at  the  interval  under  consideration 
will  contain  one  watt  of  power.  The  signal-to-noise  ratio  at 
this  interval  is  1/(1/512)  or  27.1  dB.  However,  the  signal- 
to-noise  ratio  at  the  input  to  the  processor  is  1/1  or  0  dB. 
So,  in  evaluating  the  performance  of  a  detection  method,  care 


Table  1.  Summary  of  probability  of  detection  results 
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must  be  used  in  defining  the  signal -to-noise  ratio.  All  of 
the  signal-to-noise  ratios  contained  in  this  report  are  fre¬ 
quency  interval  (bin)  signal-to-noise  ratios  based  on  a 
block  length  of  1024. 

Chapter  II  of  this  report  contains  a  detailed  discussion 
of  FABLUS  including  analysis  of  the  errors  associated  with 
application  of  the  FFT  and  the  operations  necessary  to  pro¬ 
duce  the  arrays  of  numbers  representing  the  energy  content 
of  noise  only  and  combination  of  signal  and  noise.  Chapter 
III  discusses  the  statistical  techniques  investigated  to 
determine  the  significance  of  the  spectral  components.  Experi¬ 
mental  results  obtained  using  FABLUS  are  contained  in  Chapter 
IV.  Chapter  V  presents  conclusions.  Derivations  of  analyti¬ 
cal  results  are  contained  in  Appendices  A,  B  and  C.  A  list 
of  references  is  provided. 
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II.  THE  FABLUS  PROCESS 


A.  PROCEDURE 

Since  the  FFT  processor  is  a  digital  device,  the  analog 
input  voltage  must  be  sampled  and  quantized.  Fig.  4  shows  an 
analog  waveform  x(t)  and  the  sampling  operation.  (The  errors 
caused  by  quantization  are  not  addressed  in  this  report.) 

The  sampling  interval  T  and  the  number  N  of  sample  points 
in  a  data  record  are  values  which  must  be  chosen  carefully 
if  the  FFT  processor  output  is  to  represent  the  spectrum  of 
the  input.  Appendix  A  contains  a  brief  review  of  Fourier 
transform  theory  and  the  Discrete  Fourier  Transform  (DFT) . 

The  output  of  the  FFT  processor  is  shown  in  Fig.  5  where 
X(mfQ)  represents  the  discrete  Fourier  transform  of  x(t); 
m  is  the  set  of  integers.  The  output  consists  of  N  complex 
numbers;  however,  only  the  first  N/2  are  unique.  The  spectra 
is  "folded"  about  the  N/2  value.  Therefore,  the  usable  output 
consists  of  N/2  complex  numbers  representing  the  inphase  and 
quadrature  components  of  the  transform.  These  components 
are  the  transform  values  at  discrete  frequencies  at  intervals 
of  1/NT  Hz  in  the  range  0  to  1/2T  Hz.  We  refer  to  an  interval 
of  1/NT  Hz  as  a  frequency  bin.  We  denote  the  output  numbers 
by  when  the  combination  of  signal  and  noise  is  the  input 
and  by  when  noise  only  is  the  input.  These  numbers  vary 
from  one  block  of  data  to  the  next  because  of  the  random 
nature  of  the  input  (random  phase)  relative  to  the  switching 
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Fig.  5.  Discrete  Fourier  Transform  of  sampled 
waveform. 
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and  data  truncation  which  is  part  of  the  FFT  processor  opera¬ 
tion.  To  remove  this  randomness  the  squared  magnitude  of  the 

2  2 

transform  (Re(  )  +  Im(  )  }  is  computed.  The  squared  magni¬ 

tudes  are  insensitive  to  phase  variation  and  are  proportional 
to  the  power  contained  in  the  input  time  record  at  each  fre¬ 
quency  bin.  Fig.  6  shows  the  FFT  output  magnitude  squared. 
Therefore,  the  first  processing  step  is  to  convert  the  FFT 
output  to  magnitudes  squared  (|GCI  and  |GN|  )  at  each  bin. 

Averaging  of  M  output  arrays  smooths  the  magnitude  squared 
values  and  provides  values  proportional  to  energy  for  each 
bin.  To  perform  the  averaging,  the  sum  of  the  magnitudes 
squared  must  be  computed  and  stored  after  each  block  of  input 
data  has  been  processed.  In  addition  to  storing  the  sum  of 

the  magnitudes  squared,  the  sum  of  the  squares  of  the  magni- 

2  2  2  2 

tudes  squared  ( [ | Gc |  ]  and  t  j  G^ |  ]  )  is  also  stored  to  allow 

2 

computation  of  the  variance  sc  of  the  signal  plus  noise  data 

2 

and  of  the  variance  sN  of  the  noise  only  data.  Since  the 

switching  provides  alternately  noise  only  and  a  combination 

of  signal  and  noise,  a  total  of  2M  blocks  of  data  must  be 

- 2  2 

processed  before  computation  of  means  ( | Gc |  and  |gn|  )  and 

2  2 

sample  variances  (sr  and  sM  )  can  be  done.  In  the  remainder 
of  this  report  M  will  be  referred  to  as  ensemble  size  and 
will  refer  to  the  taking  of  a  total  of  2M  samples  of  input 
voltage.  These  four  arrays  are  the  numerical  output  from 

which  the  estimate  of  the  spectrum  of  the  signal  is  formed. 

2  2  2  2 

Appendix  B  contains  derivations  of  |GC|  ,  | G^ |  ,  sc  and  s^ 

for  the  case  of  zero  mean,  white,  Gaussian  noise. 
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B.  PROCESSING  ERRORS 

The  output  of  the  FFT  processor  contains  the  information 
necessary  to  estimate  Gg(f),  the  spectrum  of  s(t).  Variations 
in  the  arrays  GN  and  Gc  due  to  random  phase  have  been  elim¬ 
inated  by  using  the  transform  magnitude  squared.  However, 
there  are  still  four  potential  sources  of  error  which  could 
cause  this  estimate  of  Gg(f)  to  differ  from  its  true  value. 
Three  of  these  error  sources  are  sampling  of  the  input  wave¬ 
form,  truncation  of  the  input  data  record  and  the  discrete 
nature  of  the  output  spectrum.  These  three  are  all  associated 
with  the  FFT  process.  The  fourth  potential  source  of  error 
is  associated  with  the  switching  action  necessary  to  produce 
the  alternating  inputs  of  noise  only  and  combination  of  signal 
and  noise. 

Appendix  A  contains  a  discussion  of  the  errors  associated 
with  estimating  a  continuous  spectrum  using  an  FFT  processor. 
As  derived  there,  the  output  of  the  FFT  processor,  with  input 
x(t)  can  be  viewed  as  composed  of  three  terms: 

X(mfQ)  =  X'(mfQ)  -  Al(mf0)  -Tr(mfQ) 

In  this  equation  X(mfg)  is  the  estimate  of  X(f),  Al(mfg) 
is  the  error  caused  by  aliasing  in  the  sampling  of  x(t), 
Tr(mfg)  is  the  error  caused  by  truncation  of  the  input  data 
and  X'(mf0)  is  the  sampled  replica  of  X(f),  the  continuous 
spectrum;  fQ  is  the  frequency  interval  =  1/NT.  If  the  input 
to  the  FFT  processor  is  noise  n(t)  only,  the  output  will  be 
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~  t mf  o ) 

When  the  input  is  signal  plus  noise,  s(t)  +  n(t),  the  pro¬ 
cessor  output  becomes: 

Gc(mf0)  =  Gc'(mfQ)  -  Alc (mf q)  -  Trc(mf0) 

The  errors  due  to  aliasing  in  the  signal  plus  noise  output 
can  be  considered  as  composed  of  components  associated  with 
the  noise  (AlN)  and  components  associated  with  the  signal 
(Alg) .  If  the  signal  and  noise  are  uncorrelated,  these  two 
components  will  be  independent  and  Alc(mfQ)  can  be  written 
as  Als(mf0)  +  Al^ (mf q) .  Errors  associated  with  the  trunca¬ 
tion  of  the  input  for  signal  plus  noise  conditions  can  also 
be  written  as  two  components  Trc(mfQ)  =  Trs(mfQ)  +  TrN(mfg) 
because  of  the  linearity  of  the  Fourier  transform.  Form 

R (mf  Q )  =  Gq  (mf q)  -  G^mfg) 

=  Gs,(mfo)  ■  Als(mf0)  ■  Trs(mfo)  +  V^V 

-  Al^mfg)  -  TrN(mf0)  -  GN'(mfQ)  +  Al^mfg) 

+  TrN(mf0) 

*  Gs'(nif0)  -  Alg  (mfg)  -  Trg(mf0) 

This  result  shows  that  the  errors  caused  by  aliasing  and 
truncation  of  the  noise  only  waveform  are  eliminated  when  the 
difference  is  formed.  Therefore,  the  computed  difference 
represents  the  spectrum  of  the  signal  as  if  it  had  been  applied 
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to  the  FFT  processor  without  the  noise.  Since  R(mf0)  con¬ 
tains  these  erroneous  frequency  components  which  are  not 
eliminated  by  the  post  transform  processing,  care  must  be 
exercized  to  minimize  their  effects.  Careful  choice  of  T, 
the  sampling  period,  and  N,  the  number  of  samples  per  block, 
will  accomplish  this. 

One  of  the  major  features  of  this  detection  scheme  is 
the  switching  of  the  antenna  into  and  out  of  the  receiver 
system  input.  This  switching  action  can  be  modeled  as  a 
squarewave  alternating  between  0  and  +1.  Call  this  waveform 
j(t).  With  s(t)  the  signal  and  n(t)  the  noise  the  FFT  pro¬ 
cessor  input  x(t)  is: 

x  ( t )  =  s  (t)  j  (t)  +  n  (t) 

The  multiplication  of  s(t)  by  j(t)  produces  discontinuities 
in  x ( t ) .  Multiplication  in  the  time  domain  corresponds  to 
convolution  in  the  frequency  domain.  Denoting  the  convolution 
operation  by  *,  we  then  have  that  the  spectrum  X(f)  of  x(t) 
is: 


X(f)  *  S(f)*J(f)  +  N  ( f ) 

The  convolution  of  J(f)  with  S(f)  produces  frequency  components 
in  X ( f )  which  are  not  in  S(f).  However,  these  spurious  com¬ 
ponents  can  be  eliminated  by  blanking  the  FFT  processor  input 
during  the  transitions  of  j(t).  This  will  require  that  the 
record  length  of  the  processor  input  be  slightly  less  them 
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the  half  period  of  j (t)  and  that  the  FFT  input  blanking  pulse 
and  j(t)  be  synchronized  to  allow  the  transients  to  decay 
sufficiently  before  the  data  is  applied  to  the  FFT  processor. 
That  is,  the  effects  of  switching  are  eliminated  by  selecting 
a  complete  block  of  data  to  be  processed  during  the  time 
interval  the  switch  is  closed. 

In  summary,  the  output  of  the  FFT  processor  is  an  array 
of  numbers  representing  an  estimate  of  the  spectrum  of  the 
processor  input.  This  estimate  is  influenced  by  the  phase 
variations  of  the  data  blocks  relative  to  the  signal,  and 
by  sampling,  truncation  and  switching.  The  random  phase 
variations  are  removed  by  computing  the  magnitude  squared. 

The  errors  associated  with  the  switching  are  removed  by 
synchronization  of  the  data  taking  and  switching  operations. 
Post  transform  processing  (differencing)  eliminates  errors 
caused  by  sampling  and  truncation  of  the  noise.  The  errors 
remaining  in  the  post  transform  processed  data  are  those 
associated  with  sampling  and  truncation  of  the  signal  s(t). 
These  are  errors  in  the  FFT  process  and  are  not  caused  by 
the  switching  action. 

Fig.  7  is  a  more  detailed  block  diagram  of  the  FABLUS 
system. 
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diagram  of  the  FRBLUS  system. 
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III.  STATISTICAL  METHODS 


The  output  of  the  FPT  processor  is  four  arrays  of  numbers 
representing  the  sum  of  the  squared  magnitudes  and  the  sum 
of  the  squared  magnitudes  squared  (fourth  power)  at  each  fre¬ 
quency  bin  for  noise  only  and  signal  plus  noise  inputs.  The 
next  step  in  the  detection  process  is  to  manipulate  these 
arrays  to  perform  a  statistical  evaluation  at  each  bin.  The 
evaluation  determines  which  bins  contain  statistically  signi¬ 
ficant  components.  Those  bins  with  significant  components 
and  the  magnitudes  of  those  components  comprise  the  spectral 
estimate  of  s(t). 

There  are  at  least  two  general  courses  of  action  available 
to  perform  the  statistical  analysis.  First,  it  can  be  assumed 
that  a  signal  is  present  and  then  evaluate  the  arrays  to 
confirm  or  deny  this  hypothesis.  Alternatively,  it  can  be 
assumed  that  no  signal  is  present  and  then  evaluate  the  arrays 
to  confirm  or  deny  this  hypothesis.  These  two  approaches  are 
not  equivalent. 

Hypothesizing  that  a  signal  is  present  requires  some 
assumption  about  the  level  of  that  signal.  This  assumed  level 
is  the  smallest  which  can  be  detected  with  the  desired  degree 
of  confidence.  Then,  the  statistical  evaluation  determines 
with  1-a  confidence  that  a  signal  of  the  assumed  level  or 
higher  exists.  This  procedure  sets  a  detection  rate  and  a 
rate  of  missed  detections. 
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Hypothesizing  that  no  signal  is  present  requires  no 
assumptions  about  the  signal  level.  The  statistical  evalua¬ 
tion  determines  with  l-o  confidence  that  no  signal  is  present. 
This  method  sets  the  rate  at  which  correct  no  signal  decisions 
are  made  and  the  false  alarm  rate.  The  relation  between  these 
hypotheses,  and  the  detection  and  false  alarm  rates  are  de¬ 
picted  in  Fig.  8. 

The  no  signal  hypothesis  requires  no  assumptions  about 
signal  level  and  setting  of  the  false  alarm  rate.  Unfortu¬ 
nately,  the  signal  detection  rate  cannot  be  determined  with¬ 
out  signal  level  assumptions. 

To  minimize  the  number  of  assumptions  required,  the  hypo¬ 
thesis  no  signal  present  is  used  in  the  signal  detection 
methods  employed  here.  This  is  in  agreement  with  general 
practice  where  control  of  false  alarm  rate  is  usually  pre¬ 
ferred  to  control  of  signal  detection  rate.  Various  methods 
are  investigated.  The  two  most  useful  methods  are  first  des¬ 
cribed  in  detail  in  the  following. 

A.  METHOD  I — DIFFERENCE  OF  TWO  MEANS 

This  method  requires  fewest  assumptions.  The  only  assump¬ 
tion  is  that  M  (the  number  of  samples  used  in  acquisition  of 

the  arrays  stored  by  the  FFT  processor)  be  large  enough  so 

2  2 

the  sample  means  (1^1  and  |GC|  )  are  approximately  Gaussian 

2 

distributed  and  the  sample  variances  (sN2  and  sc  )  are  approxi- 

2  2 

mately  the  population  variances  (aN  and  )  .  Sample  sizes 
as  small  as  10  or  even  less  can  meet  the  first  requirement 
to  reasonable  accuracy.  It  is  more  difficult  to  predict  the 
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number  of  samples  required  to  meet  the  second  condition. 
Values  of  30  or  40  are  cited  in  some  texts  [Ref.  2,3].  How¬ 
ever,  since  the  population  variances  are  unknown,  no  other 
reasonable  action  is  available  and  this  substitution  is  used 
for  all  values  of  M. 

Under  the  assumptions  listed  above,  the  random  variables 
2  2 

IGjjI  and  |GC|  can  be  considered  Gaussian  with  means  u  (the 

actual  mean  of  the  noise  population)  and  v  (the  actual  mean 

2 

of  the  signal  plus  noise  population)  and  variances  oN  /M  and 

2  2 
<?c  /M  respectively.  Form  the  random  variable  R  =  |GC| 

2 

IG-J  which  is  Gaussian  distributed  with  mean  (v-u)  and 

i  i 

2  2 

variance  ((ac  +  oN  ) /M) .  As  indicated  above  the  population 
variances  must  be  replaced  with  the  sample  variances  to  yield: 

2  2 

Random  Variable  R:  Normal  ((v-u),(sc  +sN  ) /M) 


Without  loss  of  generality  we  cam  assume  that  no  signal  is 
present  (v  =  u)  and  make  the  null  hypothesis  HO:  u  =  v. 
The  method  proceeds  as  follows: 

For  each  frequency  bin  i: 


72  2  2  2 

a.)  Compute  GN  i  ,  Gc  i  ,  sNi  and  sci 


7? 


b.)  Compute  R^  =  | Gc I i  -  | | i  and 


3Si2  =  ^SCi2  +  SNi2)/M 


c.)  Compute  Cv  »  2assi 


where  z  =  the  standardized  normal  variable  such 

a 


that  the  area  under  this  density  function  from 
z  to  »  equals  a. 


d. )  If  <  Cv  accept  HO  and  make  decision  DO  =  no 

signal  at  bin  i. 

e. )  If  Ri  >=  Cv  reject  HO  and  make  decision  Dl  =  signal 

present  at  bin  i. 

f. )  If  HO  is  accepted  set  R^  =  0;  spectral  estimate 

does  not  have  a  component  at  ith  bin. 

g. )  If  HO  is  rejected,  then  R^  is  the  estimate  of  the 

magnitude  of  the  spectral  component  of  the  signal 
at  bin  i. 

Fig.  9  is  a  block  diagram  of  this  method. 

If  no  signal  is  present  then  the  random  variable  Ri  should 

2 

have  zero  mean  with  variance  2sNi  /M  since  the  variance  is 
the  same  for  both  sets  of  M  samples.  The  probability  density 
function  for  R^  conditioned  on  no  signal  is: 

p(r|u=v)  =  /M/ (2sN/7)  exp  (-Mr^/4si;;2) 

where  r  is  a  value  of  the  random  variable  R^.  The  critical 
value  is  za/2sN//M.  The  region  to  the  left  of  the  critical 
value  represents  those  values  of  R  for  which  HO  is  accepted; 
to  the  right,  HO  is  rejected.  From  this  it  is  obvious  that 

the  choice  of  z  sets  the  false  alarm  rate. 

a 

If  a  signal  is  assumed  to  be  present,  then  v-u  =  d  and 

2  2 

R^  has  mean  d  and  variance  (sc^  +  sjji  )/M*  Tlle  probability 
density  function  for  R^  conditioned  on  a  signal  present  is: 

p(r  j  v-u  =*d)  *  /M/v/^7^) ,  exp  (-M(r-d)  2/2  (s,.2^2) ) 
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The  critical  value  is  not  the  same  as  in  the  no  signal  case 
because  the  variance  of  R  is  greater  due  to  the  signal.  The 
probability  of  detection  will  be  a  function  of  d  and  Sg. 

The  value  of  z^  is  set  by  the  false  alarm  rate.  Fig.  10  is 
a  graph  showing  both  the  no  signal  and  signal  present  density 
functions . 

To  evaluate  the  analytical  performance  of  this  spectral 
estimation  method,  we  apply  the  results  from  Appendix  B  to 
the  procedure  described  above.  The  probability  of  detection 
for  Method  I  is  P (r  > Cv)  where  Cv  is  the  critical  value. 

This  is  evaluated  in  Appendix  C  for  the  case  of  sinusoidal 
signal  and  zero  mean,  white,  Gaussian  noise  to  give: 

P  (r  >  Cv)  =  Q(za  -  (A2N2  i/M)  /W (N2A2Nn2/2  +  2N2  (n2 )  2 )  )  ) 

In  this  equation  A  is  the  peak  amplitude  of  the  sinusoidal 
~2 

signal,  n  is  the  mean  square  value  of  the  noise  process,  M 
is  the  ensemble  size,  N  is  the  number  of  points  in  the  input 
data  block  and  z^  is,  as  before,  the  value  of  the  argument  of 
the  normal  probability  density  function  such  that  the  area 
under  the  function  from  za  to  »  is  a.  The  function  Q(x)  is 
defined  as  the  area  under  the  standard  normal  density  function 
from  x  to  +®.  Fig.  11  is  a  graph  of  probability  of  detection 
versus  bin  signal-to-noise  ratio.  From  Fig.  11  we  see  that  at 
0  dB  bin  signal-to-noise  ratio,  the  probability  of  detection 
increases  by  22  percent  (from  30  to  52)  when  the  number  M  of 
records  averaged  is  doubled  (from  8  to  16) . 
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B.  METHOD  II — MEANS  ONLY 

Method  I  requires  that  the  mean  and  variance  of  both  the 

signal  plus  noise  and  the  noise  only  ensembles  be  computed. 

These  computations  can  be  time  consuming.  Method  II  was 

developed  to  reduce  the  data  and  computation  requirements. 

For  this  method  only  the  mean  values  of  the  signal  plus 

noise  and  noise  only  ensembles  are  required.  However,  the 

assumptions  required  are  more  restrictive  than  Method  I. 

In  addition  to  the  assumption  that  the  sample  means  are  Gaussian 

distributed,  it  is  assumed  that  a  linear  relationship  exists 

between  the  standard  deviation  (sVT)  of  the  noise  only  and  the 

N 

2 

mean  ( | G„ |  )  of  the  noise  only,  and  that  the  variance  of  both 

2" 

ensembles  are  equal.  That  is,  it  is  assumed  sN  =  B|Gn| 

2  2 

and  a.T  =  a „  where  B  is  a  constant. 

N  C 

It  is  difficult  to  identify  cases  where  this  linear  rela¬ 
tionship  exists  unless  analytical  methods  are  employed. 

Appendix  A  shows  that  for  Gaussian,  white  noise  these  assump¬ 
tions  are  valid  with  B  =  1.  The  equal  variance  assumption 

is  in  general  not  valid;  however,  for  small  signal  levels 
2  2 

(aS  <<:  <?N  )  >  the  error  introduced  is  small. 

72 

Under  these  assumptions ,  the  random  variables  [ |  and 
2 

|GC I  are  Gaussian  distributed  with  means  u  and  v  respectively 

2  2 

and  common  variance  a„  .  Then  the  random  variable  R .  —  I G _ I 

N  1  '  C 

— — 2 
|GN|  is  Gaussian  with  mean  (v-u)  and  variance  2c>N  /M.  With 

the  above  assumptions  we  have: 

Random  Variable  R:  Normal ( (v  -  u) ,  2B2 ( | G^ | 2 ) 2/M) 
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As  in  Method  I  we  assume  that  u  *  v  and  the  null  hypothesis 
HO:  u  *  v. 

For  each  frequency  bin  i: 

a. )  Compute  |Gc|i  and  |GN|i 

5  5 

b. )  Compute  Ri  =  |Gc|i  -  |GN|i  and 

sSi  =  ^iGjjli2/^ 

c. )  Compute  critical  value  =  Cv  =  z  s„. 

a  Si 

d. )  If  <  Cv  accept  HO  and  make  decision  DO  =  no 

signal  at  bin  i. 

e. )  If  R.  >  Cv  reject  HO  and  make  decision  D1  =  sig¬ 

nal  present  at  bin  i. 

f .  )  If  HO  is  accepted  set  R^  =  0;  spectral  estimate 

does  not  have  a  component  at  ith  bin. 

g. )  If  HO  is  rejected,  then  R^  is  the  estimate  of  the 

magnitude  of  the  spectral  component  of  the  signal 
at  bin  i. 

Fig.  12  is  a  block  diagram  of  the  method. 

Method  II  is  basically  the  same  as  Method  I.  The  only 

substantial  difference  is  the  method  of  establishing  the 

variances  and  therefore  the  critical  value. 

Assuming  no  signal  is  present  the  decision  variable  R  is 

2 

Gaussian  distributed  with  zero  mean  and  variance  2sN  /M. 
Therefore  the  probability  density  function  of  R  conditioned 
on  no  signal  is: 

P(r|u=v)  =  /M/ (28^/7)  exp(-r2M/4sN^) 
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2 

substituting  B | |  for  s 

P(r|u=v)  =  /M/(2B|GN|2/7)exp(-Mr2/(4B2{  |GNj2}2) 

The  critical  value  will  be: 

cv  *  2  /IbIg^J2/^ 

Like  Method  I  the  area  under  the  density  function  to  the  left 
of  the  critical  value  represents  the  probability  of  accepting 
HO  (correct  no  signal  decision) .  The  area  to  the  right  repre¬ 
sents  the  probability  of  incorrectly  rejecting  HO  (false  alarm) . 
If  B  is  known  exactly,  then  za  will  set  the  false  alarm  rate. 

If  a  signal  is  present  the  decision  variable  is  Gaus¬ 
sian  distributed  with  mean  (v-u) .  However,  the  variance  is 
the  same  as  in  the  no  signal  case.  Therefore,  the  density 
function  conditioned  on  signal  present  will  be: 

P (r ( v-u  =  d)  =  /M/(2B!GNj2/7)exp(-M(r-d)2/4B2{ |Gn|2}2) 

The  area  under  the  curve  to  the  left  of  the  critical  value 
represents  the  probability  of  accepting  HO  (missed  detection) . 
The  area  to  the  right  is  the  probability  of  correcting  de¬ 
tecting  the  signal.  In  this  case  the  probability  of  detection 
will  be  a  function  of  d,B  and  |gn|  .  Fig.  13  shows  the 
density  functions  under  signal  and  no  signal  conditions. 

C.  OTHER  METHODS 

In  addition  to  the  methods  described  above,  we  also  in¬ 
vestigated  two  methods  which  require  the  assumption  that  the 
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13.  Probability  density  function*  for  M*thod  II 


population  being  sampled  is  Gaussian  distributed.  The  two 
methods  are  Smith-Satterthwaite  test  for  the  difference  of 
two  means  using  Student  t  statistics  and  the  ratio  of  two 
variances  using  F  distribution  statistics.  The  results  for 
both  these  methods  are  significantly  poorer  than  those  of 
either  Method  I  or  II. 

1.  Smith-Satterthwaite  Test — Difference  of  Means 

This  test  suffered  from  two  identified  shortcomings. 
Unlike  the  two  previous  methods  where  the  accuracy  of  the 
assumptions  increased  with  increasing  M,  the  validity  of 
assuming  the  population  being  sampled  is  Gaussian  distributed 
is  not  a  function  of  M.  Therefore,  errors  caused  by  using 
this  assumption  persist  for  any  ensemble  size.  In  addition, 
because  the  critical  value  is  determined  by  calculating  the 
degrees  of  freedom  from  the  variance  values,  it  is  difficult 
to  get  consistent  false  alarm  performance  [Ref.  4] . 

2 .  F  Distribution — Ratio  of  Variances 

The  ratio  of  two  variances  test  also  is  afflicted 
with  the  errors  caused  by  assuming  the  sampled  population 
is  Gaussian  distributed.  In  addition  the  variability  of  the 
ratio  of  two  variances  is  quite  large  at  the  relatively  small 
ensemble  sizes  of  interest  here.  This  coupled  with  the 
rather  large  right  hand  tail  of  the  F  distribution  causes 
the  performance  to  be  substantially  worse  than  any  of  those 
already  discussed  [Ref.  5]. 
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D.  REFERENCE  CALCULATION 

One  of  the  serious  sources  of  error  in  Methods  I  and  II 
is  the  variability  of  the  critical  value  Cv.  As  long  as  the 
process  is  stationary  the  critical  value  should  be  a  constant: 


Cv  =  z  a  _ 
a  C 

Since  the  population  variances  are  unknown,  we  use  sc  as  an 
estimate  of  o c.  The  changing  value  of  Cv  will  cause  errors 
in  detection  in  addition  to  those  due  solely  to  the  random¬ 
ness  associated  with  the  decision  variable  R^. 

To  gain  insight  into  the  amount  of  error  caused  by  the 
changes  in  the  critical  value,  this  reference  method  was 
developed.  To  establish  a  constant  critical  value,  the  values 
of  the  variance  of  the  noise  and  signal  plus  noise  ensembles 
must  be  known  exactly.  To  obtain  a  good  estimate  of  these 
variances  the  statistics  of  both  the  noise  and  signal  plus 
noise  ensembles  were  collected  during  the  testing  of  Methods 

I  and  II.  This  information  was  then  evaluated  to  establish 

2  2 

the  values  of  ac  and  oN  based  on  all  the  data  taken  during 

2 

the  tests.  These  numbers  were  used  to  calculate  oc  .  The 
computation  proceeds  as  follows: 


At  each  frequency  bin  i: 

- j.  - , 

a . )  Compute  | Gg | i  and  | GN | ii 


b .  )  Compute  Ri  *  I Gs 1 1  “  I GN I i 

c. )  Read  in  oc  and  compute  za<?c 
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d. )  If  R^  <  z  o  then  accept  HO;  no  signal  at  bin  i 

e. )  If  R^  >  reject  HO;  signal  present  at  bin  i 

f. )  If  HO  accepted  set  R^  =  0;  spectral  estimate  has 

no  component  at  bin  i. 

Because  this  method  eliminates  the  errors  caused  by 
variability  in  the  critical  value,  the  detection  results 
should  be  the  best  obtainable  for  a  given  signal  level. 

Thus  the  detection  results  obtained  with  this  method  are 
used  as  a  reference  to  compare  the  performance  of  Methods  I 
and  II.  Fig.  14  is  a  block  diagram  of  the  computation  pro¬ 
cess  used  to  obtain  the  reference  results. 

E.  SIGNAL-TO-NOISE  RATIOS 

The  bin  signal-to-noise  ratio  is  the  ratio  of  the  average 

signal  magnitude  squared  to  the  average  noise  magnitude 
2  jr 

squared  (|Gg|  /|Gn!  )  for  a  particular  frequency  bin.  For 
the  case  considered  above  the  signal-to-noise  ratio  is,  from 
Appendix  A, 

SNR1  »  A2N2/4Nn2 

2  2 
-  A  w/4n 

This  is  a  function  of  N,  the  input  data  block  length,  when 
using  the  FFT  algorithm. 

Alternatively  an  expression  for  the  composite  signal-to- 
noise  ratio  can  be  developed  by  solving 

|  Gjj  | 2  -  Nn^  and  1GS12  *  A2N2/4 
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Fig.  14.  Block  diagram  of  the  procedure  used  to 
compute  the  reference  results. 


for  A  and  n  and  forming: 


SNR2  =  A2/2n2 

«  2|gs|2/n|gn|2 

In  this  case  the  factor  N  is  in  the  denominator  whereas  be¬ 
fore  it  was  in  the  numerator.  For  N  =  1024  the  following 
relationship  exists  between  these  alternative  representations 
when  SNR1  and  SNR2  are  in  decibels  and  where  10  log  N/2  = 

10  log  512  =  27.1. 

SNR2  =  SNR1  -  27.1  dB 

The  expression  shown  for  SNR1  is  used  throughout  this 

report  since  the  methods  employed  make  spectral  estimates  on 

a  bin-by-bin  basis.  Consequently,  the  noise  power  (or  energy) 

at  a  particular  frequency  bin  influences  only  the  spectral 

estimate  at  that  bin.  Furthermore,  in  the  general  case  of 

wideband,  noise-like  signals  and  non-Gaussian  noise,  SNR2  is 

2  ~~2 

not  calculable  because  formulae  for  A  and  n  as  functions 
of  the  average  magnitudes  squared  are  not  known. 


IV.  EXPERIMENTAL  RESULTS 


A.  NARROWBAND  RESULTS 

To  demonstrate  the  capabilities  of  this  method,  a  hard¬ 
ware  evaluation  was  conducted  using  a  signal  generator,  noise 
source,  spectrum  analyzer  and  minicomputer.  A  square  wave 
signal  was  used  to  allow  simultaneous  processing  of  multiple 
signal  levels.  The  noise  source  produced  white,  Gaussian 
noise.  The  spectrum  analyzer  did  the  sampling  and  obtained 
the  FFT.  The  minicomputer  provided  the  switching  control, 
storage  and  computational  capacity  necessary  to  carry  out  the 
statistical  analyses. 

A  total  of  1000  trials  was  conducted  with  ensemble  size 
M  of  eight.  After  each  set  of  data  was  collected  and  the 
various  sums  computed,  the  values  at  the  frequency  bins  con¬ 
taining  the  primary  and  third,  fifth,  seventh,  ninth  and 
eleventh  harmonic  frequencies  of  the  input  square  wave  were 
stored  as  well  as  data  from  five  bins  which  contained  no  sig¬ 
nal.  After  the  1000  trials  were  completed  the  data  collected 
was  analyzed  using  the  previously  described  Methods  I  and  II. 

A  preliminary  analysis  was  conducted  using  z^  =  1.96  for  Method 

I  and  z  =1.96  and  B  =  1  for  Method  II  (2.5  percent  false 
a 

alarm  rate) .  Because  of  the  variation  in  the  validity  of 
the  assumptions  used,  the  desired  false  alarm  rate  was  not 
achieved.  The  value  of  zq  was  then  adjusted  separately  for 
each  method  until  a  2.5  percent  false  alarm  rate  was  achieved. 
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Next  the  data  was  reprocessed  to  determine  the  average  values 
and  variances  at  each  bin  over  the  1000  trials.  These  aver¬ 
ages  and  variances  were  then  used  to  process  the  data  to  ob¬ 
tain  the  reference  results.  Table  2  shows  the  detection 
results  for  these  trials. 

The  data  obtained  for  the  1000  trials  with  M  =  8  was  then 
reprocessed  to  yield  data  for  500  trials  with  ensemble  size 
16.  This  data  was  then  processed  using  Methods  I,  II  and 
the  reference  procedure.  These  results  are  presented  in 
Table  3. 

The  results  for  both  the  M  =  8  and  M  *  16  trials  are 
plotted  in  Fig.  15  along  with  the  analytical  probability 
of  detection  curves. 

A  comparison  of  the  analytical  and  experimental  results 
reveals  that; 

1.  The  analytical  and  experimental  results  for 
Method  I  agree  very  closely. 

2.  The  reference  results  are  substantially  better 
than  either  Method  I  or  II. 

3.  Method  I  is  slightly  better  than  Method  II. 

The  very  close  agreement  between  the  analytical  and  ex¬ 
perimental  results  for  Method  I  indicates  that  the  assumption 
used  (sample  mean  is  Gaussian  distributed)  is  generally  valid 
for  the  case  of  white,  Gaussian  noise  and  sinusoidal  signals. 
The  results  are  not  expected  to  change  for  other  signal  for¬ 
mats,  since  any  waveform  encountered  in  practice  can  be 
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TABLE  2.  Probability  of  detection  Ensemble  size  3 

false  alarm  rate  2.5';. 


PROBABILITY  OF  DETECTION  <X> 

BIN 

SNR  METHOD  I  METHOD  II  REFERENCE 

12.2  100  100  100 

2.6  54.2  50.3  32.4 

-2.2  13.2  11.7  23.9 

-4.8  6.7  6.3  11.6 

-7.2  3.3  4.1  6.3 


-9.7 


2.8 


4.4 


4.6 


TABLE  3.  Probability  of  detection  <50.  Ensemble  size  16, 
false  alarm  rate  2.55;. 


PROBABILITY  OF  DETECTION  <50 


BIN 

SNR 

METHOD  I 

METHOD  II 

1 

REFERENCE 

12.2 

100 

108 

100 

2. 6 

83.2 

34.6 

99.3 

-2.2 

29.6 

22.2 

56.0 

-4.3 

12.6 

12.6 

23.6 

-7.2 

4.3 

4.2 

11.2 

1 

-9.7 


5.4 


6. 2 


7.0 


expressed  as  a  sum  of  sinusoids.  The  validity  of  the  assump¬ 
tion  for  other  types  of  noise  is  open  to  question.  However, 
the  ability  of  the  method  to  detect  the  signals  will  not  be 
greatly  changed  by  the  presence  of  non-Gaussian  noise  if  the 
ensemble  size  is  large  enough  to  allow  the  Gaussian  dis¬ 
tributed  sample  mean  assumption. 

As  expected  the  reference  results  show  greater  detection 
ability  at  all  signal  levels.  They  are  included  here  only 
to  indicate  an  upper  limit  of  detection  under  these  conditions. 
It  may  be  difficult  to  implement  such  a  method  since  it  re¬ 
quires  prior  knowledge  of  the  statistics  of  the  processes 
involved. 

A  comparison  of  the  probabilistic  characteristics  of 

Methods  I  and  II  (see  Figs.  10  and  13)  indicates  that  Method 

II  should  have  higher  probability  of  detection  than  Method  I 

because  the  critical  value  will  be  less  under  otherwise 

identical  conditions.  However,  the  experimental  results  show 

precisely  the  opposite.  The  smaller  critical  value  obtained 

using  Method  II  is  due  to  the  use  of  the  noise  only  data  to 

compute  the  variance  of  R.  The  resulting  smaller  critical 

value  does  increase  the  probability  of  detection,  but  it  also 

increases  the  false  alarm  rate.  When  the  value  of  z  is  ad- 

a 

justed  to  regain  the  desired  false  alarm  rate,  the  probability 
of  detection  decreases  to  the  levels  shown  in  Tables  2  and  3. 
The  increase  in  the  false  alarm  rate  is  not  predicted  by  the 
analytical  analysis  but  is  due  to  additional  noise  entering 
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the  experimental  system  along  with  the  signal.  This  type 
of  error  can  be  encountered  in  any  method  where  the  critical 
value  is  calculated  from  noise  only  data.  The  use  of  both 
signal  plus  noise  and  noise  only  data  to  calculate  the  criti¬ 
cal  value,  as  in  Method  I,  provides  a  limited  amount  of 
protection  against  this  type  of  error. 

The  marginal  increase  in  detection  capability  which 
Method  I  shows  over  Method  II  is  probably  due  to  the  greater 
amount  of  information  used  in  the  Method  I  process.  Since 
Method  I  uses  both  the  means  and  variances  to  make  spectral 
estimates,  the  extremes  in  the  means  may  be  counterbalanced 
to  some  extent  by  the  variances.  This  might  produce  a  smaller 
overall  fluctuation  in  the  decision  variables.  Method  II 
uses  the  means  to  set  the  critical  value  and  to  compute  the 
value  of  the  decision  variable.  Therefore,  this  process 
would  be  more  sensitive  to  the  variations  of  the  means. 

However,  the  slightly  increased  detection  capability  of 
Method  I  does  not  occur  without  cost.  Implementation  of 
Method  I  requires  additional  computation  time  since  the  arrays 
of  squared  magnitudes  squared  must  be  computed  and  requires 
twice  the  storage  capacity  since  four  arrays  must  be  stored 
instead  of  two. 

A  graph  of  ensemble  size  versus  signal -to-noise  ratio  for 
50  and  95  percent  probability  of  detection  is  included  as 
Fig.  16.  Additionally,  Fig.  17  is  a  graph  of  probability  of 
detection  versus  bin  signal-to-noise  ratio  for  2.5,  5  and  10 
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IL1TY  OF  DETECTION 


percent  false  alarm  rates.  Figs.  16  and  17,  which  are  ob¬ 
tained  from  the  equations  derived  in  Appendix  B,  indicate 
the  trade-offs  involved  using  these  detection  methods. 

B.  WIDEBAND  RESULTS 

The  FABLUS  process  makes  decisions  on  a  bin-by-bin  basis, 
and  so  it  can  be  applied  to  wideband  signals  having  continu¬ 
ous  spectra.  To  test  this  capability,  a  pseudo-random  se¬ 
quence  of  length  1,048,575  was  used  as  a  signal.  Because  a 
noise-like  signal  is  used,  the  signal  spectrum  can  only  be 
specified  on  an  average  or  mean  basis.  Fig.  18a  shows  the 
spectrum  of  the  pseudo-random  sequence  based  on  an  average 
over  500  data  blocks.  Fig.  18b  shows  the  spectrum  of  the 
noise  also  based  on  500  samples.  The  overall  signal-to- 
noise  ratio  is  0  dB.  The  bin  signal-to-noise  ratio  has  maxi¬ 
mum  values  of  7.5  dB  in  the  main  lobe  and  -6.8  dB  in  the 
first  sidelobe. 

For  a  given  ensemble  size  the  spectral  estimate  with  noise 
can  be  no  better  than  the  signal  only  spectral  estimate  based 
on  the  same  number  of  samples.  Fig.  19a  shows  the  average 
of  the  signal  only  spectrum  when  M  is  8,  and  Fig.  19b  shows 
the  spectral  estimate  when  noise  ij  added  for  ensemble  size 
M  *  8.  Figures  20a  and  20b  are  these  results  when  M  is  16. 

As  expected  the  estimate  for  M  =  16  is  significantly  bet¬ 
ter  than  for  M  ■  8.  In  addition  to  showing  less  gaps  in  the 
main  lobe  area,  the  M  *  16  estimate  has  considerably  less 
variation  in  the  amplitudes  of  the  spectral  components. 
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20.  Measured  and  estimated  spectra;  ensemble  size  ot  16.  composite 


Neither  estixnate  shows  significant  components  in  the  side- 
lobe  area.  This  is  caused  by  the  low  bin  signal-to-noise 
ratio  in  the  sidelobe  area.  At  the  peak  of  the  sidelobe  the 
bin  signal-to-noise  ratio  is  -6.8  dB.  At  this  value,  the 
probability  of  detection  is  7.5  percent  for  M  =  16. 
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V.  CONCLUSIONS 


The  FABLUS  process  provides  an  estimate  of  the  spectrum 
of  an  unknown  signal.  The  process  performs  equally  well  with 
narrowband  and  wideband  signals.  Two  implementations  of  the 
FABLUS  process  are  evaluated.  One  method  requires  twice  the 
memory  capacity  and  additional  processing  time,  but  it  pro¬ 
duces  marginally  better  probability  of  detection,  and  it 
requires  the  least  restrictive  assumptions  about  the  signal 
spectrum.  The  accuracy  of  the  estimate  of  the  signal  spec¬ 
trum  is  limited  only  by  the  ensemble  size  and  the  inaccuracies 
inherent  in  the  FFT  process. 


APPENDIX  A 


THE  DISCRETE  FOURIER  TRANSFORM 

1 .  Equations 

A  basic  Fourier  analysis  result  is  that  a  periodic  time 
function  x(t)  can  be  represented  by  a  sum  of  sinusoids  of 
proper  frequencies  and  amplitudes.  This  sum  is  the  Fourier 
series  [Ref.  6] . 

x(t)  =  Aq  +  l  (AjjCOS  (2nnt/T)  +  BNsin  (2imt/T)  ) 
n 

where  T  is  the  period  of  x(t) .  If  x(t)  is  not  periodic  (or 
T  ->  <»)  ,  then  its  frequency  description  becomes  the  Fourier 
transform  represented  mathematically  as  [Ref.  7] : 

oo 

F  (x ( t) )  =  X(f)  =  /  x(t)  e~^2irftdt 

—  oo 

The  discrete  Fourier  transform  (DFT)  uses  sample  values 
of  x(t)  obtained  over  a  finite  length  of  time  (record  length) 
The  bandlimited,  continuous  time  function  x(t)  having  Fourier 
transform  X(f)  is  first  sampled  at  intervals  of  T  seconds  to 
produce  the  sampled  time  function  x(mT)  where  m  is  the  set 
of  integers.  Now,  let  a  finite  record  of  a  sampled  time  func 
tion  be  x(iT) ,  i  *  0,...,N-1.  The  Discrete  Fourier  transform 
of  x(iT)  is 

X (q)  *  Jx(iT) [cos (ia)  -jsin(ia)] 

i 
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where  a  =  2irq/N,  j  =  /^l,  q  =  0,...,N-1.  Here,  i  is  a 
subset  of  m,  the  set  of  integers.  Dropping  the  T  for  con¬ 
venience  and  rearranging  gives 

X (q)  =  l [x (i) cos (ia)  -  jx (i) sin (ia) ] 

i 

which  has  real  and  imaginary  parts 

Re[X(q)]  =  £x(i)cos(ia) 
i 

Im[X(q)]  =  £x(i) sin(ia) . 

i 

The  squared  magnitude  of  the  transform  which  is  used  in  this 
report  is 

1  X  (eg)  |  2  =  [ Jx(i) cos (ia)  ]  2  +  [£x  (i)  sin  (ia)  ] 2 

i  i 

=  JJx(i) x( j)cos (ia) cos ( ja)  + 

ij 

££x (i) x( j) sin (ia) sin( ja) 

ij 

=  ££x(i) x( j) [cos (ia) cos ( ja)  + 

ij 

sin (ia) sin ( ja) ] 

2 .  Errors 

From  the  sampling  theorem  [Ref.  8],  x(mT)  contains  all  the 
information  necessary  to  reconstruct  x(t)  if  T  is  selected 
so  that  1/T  is  at  least  twice  the  highest  non-zero  frequency 
component  of  X(f).  If  X(f)  exists  for  f  >  1/2T,  aliasing 


occurs  which  means  the  frequency  components  of  x(mT)  are  not 
the  same  as  those  of  x(t) .  This  is  the  first  source  of  error 
in  the  DFT . 

Since  x (mT)  contains  an  infinite  number  of  values,  for 
computational  purposes  we  must  truncate  x (mT)  to  some  reason¬ 
able  length.  This  truncation  is  accomplished  by  selecting 
some  number  N  of  samples  to  be  considered.  Therefore,  the 
record  length  Tq  is  NT.  The  truncation  process  is  also 
called  windowing  because  it  is  equivalent  to  multiplying  x(mT) 
by  a  window  function.  In  the  frequency  domain  the  multipli¬ 
cation  of  these  two  time  functions  is  equivalent  to  convolu¬ 
tion  of  the  spectrum  of  x(mT)  with  the  spectrum  of  the  window 
function.  Therefore,  the  spectrum  of  x(mT)  is  altered  by  the 
window  function.  This  is  the  second  source  of  error  in  the 
DFT  process. 

To  obtain  the  spectrum  of  x(mT) ,  the  DFT  assumes  the  avail 
able  record  repeats  and  calculates  the  Fourier  transform 
using  the  sample  values.  This  results  in  values  for  N  terms 
at  frequencies  which  are  multiples  of  1/TQ  *  1/NT  =  fQ.  This 
discrete  spectrum  X(rafQ)  may  mask  the  maximum  and  minimum 
values  of  X(f),  the  true  spectrum  of  x(t). 

The  DFT  output,  then,  consists  of  N  complex  numbers  repre¬ 
senting  the  magnitudes  and  phases  of  the  transform  at  the  N 
frequencies  spaced  1/NT  Hz  apart.  However,  only  N/2  values 
of  the  transform  are  unique.  The  values  from  N/2  to  N  repre¬ 
sent  negative  frequency  components  [Ref.  9).  Viewed  from 
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another  perspective,  it  is  unreasonable  to  expect  a  linear 
process  like  the  DFT  to  produce  2N  linearly  independent  out¬ 
puts  (real  and  imaginary  parts  at  N  frequencies)  with  only 
N  linearly  independent  inputs . 

We  can  summarize  these  various  effects  by  comparing  the 
DFT  output  with  a  sampled  version  of  X(f) . 

Let  X' (mfg)  =  the  sampled  version  of  X(f) 
and  let  X(mf^)  =  the  DFT  output. 

Then,  X(mfQ)  *  X' (mfQ)  -  Al(mfQ)  -  Tr(mfQ)  where  Al(mfQ) 
represents  the  distortion  in  Xfinf^)  caused  by  aliasing  and 
Tr(mfg)  represents  the  distortion  in  X(mf0)  caused  by  trunca 
tion.  The  aliasing  and  truncation  effects  can  cause  X(mfQ) 
to  increase  or  decrease. 

The  FFT  is  a  computationally  efficient  algorithm  used  to 
generate  the  DFT  values. 
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APPENDIX  B 


DERIVATION  OF  STATISTICAL  PROPERTIES  OF  SINUSOIDAL  SIGNAL 
AND  UNCORRELATED,  GAUSSIAN  NOISE 

In  this  appendix,  we  derive  the  equations  for  the  mean 
and  variance  of  the  frequency  bin  squared  magnitudes  which 
are  computed  by  the  FABLUS  process.  The  derivatic-  considers 
the  cases  of  noise  only  and  sinusoid  plus  noise.  That  is,  we 
derive  here  E  £ | GN (p) | 2 ] ,  Var [ | (p)  | 2 ] ,  E [ | GQ (p)  | 2 ]  and 
Var[|Gc(p)  |2]  . 

1.  Noise  Only 

From  Appendix  A,  if  we  let  x(i)  =  n(i)  =  noise  only  and 
let  Gn  represent  the  transform  under  noise  only  conditions, 
then 

2 

1  Gn  ( p )  [  =  ££n(i)n(j)  [cos  (ia)cos  ( ja)  +  sin(ia)  sin(  ja)  ] 

The  expected  value  of  the  squared  magnitude  becomes: 

E{  |  Gj.  (p)  j  }  =  ££n(i)n(j)  [cos  (ia)  cos  ( ja) +sin(ia)  sin  ( ja)  ] 

ij 

If  n(i)  is  independent  from  sample  to  sample  and  recognizing 
that  the  only  random  quantities  are  the  n(i),  then  the  ex¬ 
pected  value  is  really  the  result  of  an  ensemble  average  which 
becomes  (assuming  stationary  noise) . 

n(i)n  ( j)  *  0  for  i  ?  j 
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E{|GN(p)|2}  =  £n  (i)  n  (I) [cos2 (ia)  +  sin2(ia)J 

i 

=  Nn (i) 2 
~ 

=  Nn 

Now,  the  squared  value  of  the  magnitude  squared  is 
[ j  X (p)  [ 2] 2  =  { Ilx(i) x( j)  [cos (ia)cos ( ja) +sin(ia) sin( ja) ] }2 

ij 

=  x(j)  x(k)x(l)  [cos (ia) cos ( ja) + 

ijkl 

sin (ia) sin ( ja) ] [cos (ka) cos (la)  + 

sin(ka) sin(la) ] 

If  x(i)  =  n(i)  =*  noise  only,  then  the  mean  square  value 
becomes : 


If  n(i)  is  independent  from  sample  to  sample  and  again  recog¬ 
nizing  that  the  n(i)  are  the  only  random  quantities,  then 
this  expression  has  values  only  when  i*j»k*lori*j 
and  k  =  1  or  i  =  k  and  j  *  1  or  i  *  1  and  j  *  k. 
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There  will  be  N  conditions  when  i  =  j  =  k  =  1. 


E { [ i Gn ( p ) ( 2 3  2 }  =  £n(i) 4 [cos2 (ia) +sin2 (ia) ]  [cos2 (ia) +sin2 (ia) ] 


+  ££n(j)n{ jjn(k)n(k) [cos2 ( ja) +sin2 ( ja) ] [cos2 (ka) +sin2 (ka) ] 
jk  j?*k 

2 

+  JJn(  j)n(k)n(  j)n(k) [cos ( ja) cos (ka) +sin ( ja) sin (ka) ] 
jk  j*k 

^ _ _ _ r-rr-  2 

+  ££n ( j) n (k) n (k)n  ( j)  [cos  ( ja)  cos  (ka)  +sin ( ja)  sin  (ka)  ] 
jk 


The  cosine  squared  plus  sine  squared  terms  are  unity  and  the 
use  of  the  produce  trionometeric  identity  yields: 


E  {  [  |  G  (p)  |  2  3  2  =  £n  (i)  4  +  nn(j)n(j)n(k)n(k) 

N  i  jk  j*k 

+  [n  (j)n(k)n  ( jjn  (k)  ]  +  (n(  j)n(k)n(k)n(  j)  ] }  (cos2  [  ( j-k)  a]  } 
jk  j^k 

-  Jn (i) 4  +  J£n(j)2n(k)2  + 

i  jk 

2ZIn(j)2n(k)2cos2[(j-k)a] 
jk  j/k 


E(  [  |G„(p)  |  2]2}  =  Nn4  +  (N2  -N)  (n2)2  +  2(n2)2ncos2[(j-k)a] 

N  jk 

=  Nn4  +  (2N2  -3N)  (n2)2 


Now,  Var(x)  *  E(x2)  -  [E(x)]2.  Therefore, 


Var(|GN(p)  |2]  -  Nn*  +  (2N2  -  3N)  (n^)  2  -N2(r?)2 

~T  2  ~ 7  2 

-  Nn  +  (N*  -  3N)  (n  ) 
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4  2  2 

If  n(i)  is  Gaussian  distributed  then  n  =  3 (n  )  and  so 
2  2  J  2 

VartlG^p)  |  ]  *  N  (n  )  which  is  a  result  used  to  define  a 
signal-to-noise  ratio  in  Chapter  III. 

2.  Signal  Plus  Noise 

If  x(i)  =  s(i)  +  n(i)  =  Signal  plus  Noise,  then  letting 
Gc  represent  the  transform  under  signal  plus  noise  conditions 
the  expected  value  of  the  magnitude  squared  becomes: 

E [ | Gc (p) t  2 ]  *  ll  [s  (i)+n  (i) ]  [s ( j) +n ( j) ] (cos (ia)cos (ja) 

ij 

+sin (ia)sin (i j) T 

«  ££ {s (i) s ( j ) +s  (i) n ( j ) +s ( j )  n  (i)  +n  (i)  n ( j ) }  [cos”(ia) 

ij 

cos ( ja) +sin (la) sin (i j ) J 


If  s(i)  and  n(i)  are  independent,  n (ij  =*  0  and  recognizing 
that  the  sine  and  cosine  terms  are  not  random  gives 


E [  | G_ (p)  |  2 ]  =  ^sTUiTITlcos  (ia)  cos(ja)+sin(ia)sin(  ja)  ] 

C  ij 

+  ^n  (i)  n  ( j  )  [cos  (ia)  cos  ( ja)  +sin  (ia)  sin  ( ja)  J 
ij 

— T~~ 

From  previous  work,  the  second  term  equals  Nn  (i) .  If 

s(i)  =  Acos  (2irfci/N) ,  then  s(i)  is  deterministic,  and  we  have 

by  letting  F  =  2irfc/N 


E[|Gc(p)|2] 


^cos  (Fi)  cos  (Fj)  [cos  (ia)  cos  ( ja)  +sin(ia)  sin(  ja)  ] 

ij 


+  Nn2 
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When  F  f  p  this  reduces  to  the  noise  only  case. 


For  F  =  p  we  have 


£ [ j  Gr (p)  j  2  3  =  A2N2/4  +  Nn2 


Now  consider  the  value  of  the  squared  magnitude  squared: 


I  |X(p)  2|]  2  =  IHIx(i)x(  j)  x(k)  x(l)  [cos  (ia)  cos ( ja)+sin(ia)  sin  ( ja)  ] 
ijkl 

[cos (ka) cos (la) +sin (ka) sin (la) ] 


For  x(i)  =  s(i)  +  n(i)  we  get: 


2,2 


[i6r(P)n 


^^[3(i)+n(i)  ]  [ s ( j ) +n  ( j )  ]  [s(k)+n(k)  ]  [s(l)+n(l)  ] 


[cos (ia) cos ( ja) +sin (ia) sin ( ja) ] 


[cos  (ka)  cos  (la)  -tsin(ka)  sin  (la)  ] 


Now,  the  expected  value  of  the  squared  magnitude  squared  becomes : 


E([|Gp(p)  |  2]  2 


imCs(i)s(j)s  (k)  s  (1)  +  s(i)s(k)n(j)n(l)  + 

ijkl  _  _ 

n  (i)  s  ( j )  s  (k)  s  (1)  +  n(i)n(  j)s(k)s  (1)  + 


s  (i) s(  j)n(k)  s  (1)  +  s(i)n(j)n(k)s(l)  + 
n(i)s(  j)n(k)s(l)  +  n  (i)  n  ( j)  n(k)  s  (1)  + 


s (i) s ( j ) s (k)n  (1)  +  s (i) s (k)n ( j)n(l)  + 
n  (i)  a  ( j)  s  (k)n (if  +  n(i)n(  j)  s(k)n(l)  + 
s  (i)  s  ( j)n(k)nd)  +  s  (i) n'(j )  n(k)  s (l)  + 
n(i)s(j)n(k)n(l)  +  n(i)n(  j)n(k)n(D  1 
[cos (ia) cos ( ja)  +  sin (ia) sin ( ja) T 
[cos (ka) cos (la)  +  sin (ka) sin (la) ] 
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If  n(i)  has  zero  mean  and  zero  third  moment,  and  if  n(i)  is 
independent  of  s(i),  and  if  s(i)  is  deterministic,  and  recog¬ 
nizing  that  only  n(i)  is  random,  then 

E{  [  |Gr(p)  I  2]  2}  =  mi  [s (i) s ( j) s  (k) s  (1)  +  n(iTn  (5)s(kis(l)  + 

u  ijkl 

s(i)n(j)n(k)s(l)  +  n  (ITs  ( j )  n  (k]  s  (l)  + 
s  (i)  s  (k)n(  j)n{l)  +  n  (iYs  ( j)  s  (k)n(l)  + 
sTi)s(j)n(k)n(lj  +  n(i7 n (j) n ( k ) n (l )  ] 
(cos (ia) cos ( ja)  +  sin  (ia) sin ( ja) ] 

(cos (ka) cos (la)  +  sin(ka) sin(la) ] 

Assuming  n(i)n(j)  =  0  for  i  ?  j  (sample  to  sample  independence), 
then 


E{  [  J G_ (p)  |2]2}  =  III[s(i)s(j)s(k)s(l)  (cos(ia)cos(ja)  + 

ijkl 

sin (ia) sin ( ja) ] [cos (ka) cos (la)  + 
sin (ka) sin (la) ] 

+  2£££n(i).  s ( j) s (k) (cos2 (ia)  +  sin2(ia)J 
ijk 

[cos ( ja) cos (ka) +sin ( ja) sin (ka) ] 

+  4[I[n(i)  s  (j) s  (k)  [cos  (ia) cos  (ja)  + 
ijk 

sin(ia)sin(  ja)  (cos(ka)cos  (ia)  + 
sin(ka)  sin(ia)  ] 

+  ^JJn(i)n(j)n(k)n(l)  [cos (ia) cos ( ja)  + 
sin(ia)sin( ja) ] (cos(ka)cos(la)  + 
sin(ka) sin (la) ] 
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Assuming  s(i)  =  Acos(Fi)  where  F  =  2 tt f c/N  and  n(i)  is  Gaussian, 
then 

E  UGc(p)|2]2  «  A4N4/16  +  2[A2N3!?/4]  +  4[A2N3n2/8] 

2  T  2 
+  2n (n  j * 

4  4  2  3” 2"  2  ~2  2 

-  A^N^/16  +  A  n  n  /2  +  (n^) z 

Now,  Var(x)  *  E[x2l  -  [E(x)J2.  Therefore, 

Var[|Gc(p) 1212  =  A2N3n^/2  +  N2(n2)2 
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tv  .  "V 


V 


APPENDIX  C 


PROBABILITY  OF  DETECTION  FOR  SINUSOIDAL  SIGNAL  AND 
UN CORRELATED,  GAUSSIAN  NOISE 

Using  the  assumptions  for  Method  I  in  Chapter  III,  the 
probability  of  detection  is 

Pd  =  P  (r  >  Cv) 

where  r  is  the  decision  variable  and  Cv  is  the  critical  value. 
From  Appendix  B  for  the  case  of  white,  Gaussian  noise  and 
sinusoidal  signal, 

|Gn  !?  =  Nn^, 

Var ( | GN | 2 )  =  N2(n2)2, 

|G(;|i  =  A2N2/4  +  N?, 

Var(jGc|2)  =*  A2N3n2/2  +  N2  ( ^)2 . 

Therefore,  the  decision  variable  become 

R  -  -  |GN|2 

2  2  ~ y  ~T 

=  A  N  /4  +  Nn  -  Nn 

=  A2N2/4 

and 
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(A  «  n  /2  +  Nz(n  )*  +  N  (a  )‘)/M 


(A2N3?  +  4N2(n2)2)/2M 


=  y ( (A2N3n 


2  +  4N2(r?)2)/2M) 


The  critical  value  is: 


Cv  =  z  sc 
a  S 


=  Z 


a’/(wVn 


2  +  4N2(n2)2)/2M) 


Using  the  above  expressions, 


2%,2 


Pd  =  P(rjv-u=AN  /4) 


=  Q((Cv  -  R)/Sg} 

^2N2  /  (4  V(N3A2?/2  +  2N2(?)2))) 


=  Q ( “  (A 


y  -ft- 
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